using Plots
using Graphs
using Random
using GraphPlot
using StatsBase
using DifferentialEquations
using Base.Threads

Lucas Schmidt Ferreira de Araujo

function watts_strogatz_2d(N::Int, k::Int, β::Float64)
    L = Int(sqrt(N))  # Assuming N is a perfect square
    G = SimpleGraph(N)
    function node_index(i, j)
        return ((i - 1) % L) * L + ((j - 1) % L) + 1
    end
    for i in 1:L
        for j in 1:L
            v = node_index(i, j)
            neighbors = [
                node_index(i + di, j + dj)
                for (di, dj) in [(-1,0), (1,0), (0,-1), (0,1)]  # 4 nearest neighbors
            ]
            for u in neighbors
                add_edge!(G, v, u)
            end
        end
    end
    additional_k = (k - 4) ÷ 2
    for i in 1:L
        for j in 1:L
            v = node_index(i, j)
            for d in 1:additional_k
                add_edge!(G, v, node_index(i, j + d))  # Horizontal
                add_edge!(G, v, node_index(i + d, j))  # Vertical
            end
        end
    end
    for edge in edges(G)
        if rand() < β
            u, v = src(edge), dst(edge)
            rem_edge!(G, u, v)
            new_v = rand(setdiff(1:N, neighbors(G, u)))  # Choose new node
            add_edge!(G, u, new_v)
        end
    end
    return G
end

watts_strogatz_2d (generic function with 1 method)

function BassDiffusion(nit::Int , states0::Vector{Int} , N::Int , k::Int , β::Float64 , p::Float64 , q ::Float64)
    G = watts_strogatz_2d( N, k, β)
    states = zeros( nit , N )
    states[1,:] = states0
    for t = 1:nit-1
        states[t+1 , : ] = states[t , :]
        for j = 1:N
            n = sample( 1:N )
            state = states[t,n]
            nb = neighbors( G , n )
            if( state == 0 && length(nb) != 0)
                states_ = states[ t , nb ]
                F = sum(states_) / length( nb )
                pf = p + q*F
                states[t+1 , n] =  sample( [0,1] , Weights([1-pf,pf]) )
            end
        end
    end
    return states
end

function Adoption_State(nit::Int , states0::Vector{Int} , N::Int , k::Int , β::Float64 , p::Float64 , q ::Float64 , nsim::Int = 500)
    simulations = [ sum(BassDiffusion( nit , states0 , N , k , β , p , q),dims = 2)[:] for sim in 1:nsim ]
    return mean( simulations , dims = 1)[1]
end

function Time_Stabilization( nit::Int , states0::Vector{Int} , N::Int , k::Int , β::Float64 , p::Float64 , q ::Float64 , nsim::Int = 500 )
    time = Threads.Atomic{Float64}(0.0)
    @Threads.threads for sim = 1:nsim
        simulation = sum(BassDiffusion( nit , states0 , N , k , β , p , q),dims = 2)[:]
        if( simulation[end] == 100 )
            atomic_add!( time , findfirst( x-> x==100.0 , simulation) / nsim ) 
        end
    end
    return time[]
end

Time_Stabilization (generic function with 2 methods)

N = 100
k = 8 
nit = 300
β = range(0,1,4)
p = range(.001,.01,4)
q = range( .25 , .5 , 4);

Case I

states0 = zeros( Int, N );

plots = [ plot() for q_ in q ] 
β_ = β[1]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800) )

plots = [ plot() for q_ in q ] 
β_ = β[2]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β_ = β[3]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β_ = β[4]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β = range(0,.7,100)
nit = 500

for i in eachindex( q ) 
    plots[i] = plot()
    for j in eachindex( p )
        time = Time_Stabilization.(nit , Ref(states0) , N , k ,β, p[j] , q[i] )  
        plots[i] = plot!( β,time , label =  "p = $(p[j])")
        title!("q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800) , suptitle = "Stabilization Time")

Case II

nit = 300
states0 = zeros( Int, N )
states0[rand(1:N,8)] .= 1;

plots = [ plot() for q_ in q ] 
β_ = β[1]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β_ = β[2]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β_ = β[3]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β_ = β[4]

for i in eachindex( plots ) 
    plots[i] = plot()
    for j in eachindex( p )
        plots[i] = scatter!( Adoption_State(nit , states0 , N , k ,β_, p[j] , q[i]) , label =  "p = $(p[j])" , alpha = .1)
        title!("β=$(β_) , q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800))

plots = [ plot() for q_ in q ] 
β = range(0,.7,100)
nit = 500

for i in eachindex( q ) 
    plots[i] = plot()
    for j in eachindex( p )
        time = Time_Stabilization.(nit , Ref(states0) , N , k ,β, p[j] , q[i] ) 
        plots[i] = plot!( β,time , label =  "p = $(p[j])")
        title!("q=$(q[i])")
    end
end
plot( plots... , layout = (2,2) , size = (1600,800) , suptitle = "Stabilization Time")

• Part 3

function ODE!(du, u, P, t)
    F = u[1]
    p, q = P 
    du[1] = (1 - F) * (p + q * F)
end

function solve_model( p::Float64 , q::Float64 )
    tspan = (0.0,40.0)
    u0 = [0.0]
    ode = ODEProblem(ODE!, u0, tspan, [p,q])
    sol = solve(ode,RK4(),saveat=1.0)

    t = [ solution[1] for solution in sol.t ]
    y = [ solution[1] for solution in sol.u ];
    return t , y
end

solve_model (generic function with 1 method)

plots = [ plot() for q_ in q ] 

for i in eachindex( q ) 
    plots[i] = plot()
    for j in eachindex( p )
        t , y = solve_model( q[i] , p[j] )
        plots[i] = plot!( t , y , label =  "p = $(p[j])")
        title!("q=$(q[i])")
    end
end
plot( plots...  , size = (1600,800) )